The main purpose for this assignment’s report is to find out the spatial distribution of COVID-19 new cases in London borough and to quantify the effect of ‘Lockdown’ policy on the suppress of the increase of new cases based on statistical analysis and spatial autocorrelation.
#import library
library(sf)
library(tmap)
library(tmaptools)
library(tidyverse)
library(here)
library(janitor)
library(RColorBrewer)
library(sp)
library(spdep)
After completing loading libraries, loading all datasets needed for the analysis:
#read all dataset for the analysis
#London borough shapefile
Londonborough <-
st_read(
here::here(
'data',
'statistical-gis-boundaries-london',
'ESRI',
'London_Borough_Excluding_MHW.shp'
)
) %>%
st_transform(., 27700)
## Reading layer `London_Borough_Excluding_MHW' from data source `/Users/ziyichen/Desktop/SDSV/GIS/Data_for_CW/Code&data/data/statistical-gis-boundaries-london/ESRI/London_Borough_Excluding_MHW.shp' using driver `ESRI Shapefile'
## Simple feature collection with 33 features and 7 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: 503568.2 ymin: 155850.8 xmax: 561957.5 ymax: 200933.9
## projected CRS: OSGB 1936 / British National Grid
Read csv file–>cumulative new cases in each time period, raw csv files are in Github repo
#cumulative new cases in 3.23-4.22 period
covid_323_422 <-
read_csv(
'/Users/ziyichen/Desktop/SDSV/GIS/Data_for_CW/Code&data/data/new_cases-323-422.csv',
na = c("NA", "n/a")
) %>%
clean_names()
##
## ─ Column specification ─────────────────────────────────────────────────
## cols(
## area_code = col_character(),
## new_cases = col_double(),
## `population(10k)` = col_double(),
## `new_cases per 10k population` = col_double()
## )
#cumulative new cases in 5.23-6.22 period
covid_523_622 <-
read_csv(
'/Users/ziyichen/Desktop/SDSV/GIS/Data_for_CW/Code&data/data/new_cases-523-622.csv',
na = c("NA", "n/a")
) %>%
clean_names()
##
## ─ Column specification ─────────────────────────────────────────────────
## cols(
## area_code = col_character(),
## new_cases = col_double(),
## population_10k = col_double(),
## `new_cases_per_10k_523-622` = col_double()
## )
names(covid_323_422)
## [1] "area_code" "new_cases"
## [3] "population_10k" "new_cases_per_10k_population"
names(covid_523_622)
## [1] "area_code" "new_cases"
## [3] "population_10k" "new_cases_per_10k_523_622"
Connect the CSV file data and shapefile
#join data with shapefile and new cases data in 3.23-4.22 period
Londonborough_merged <-
Londonborough %>% left_join(covid_323_422, by = c('GSS_CODE' = 'area_code')) %>%
distinct(GSS_CODE, NAME, new_cases,new_cases_per_10k_population) #choose which column to be used in the following analysis
#join data with shapefile and new cases data in 5.23-6.22 period
Londonborough_merged_523 <-
Londonborough %>% left_join(covid_523_622, by = c('GSS_CODE' = 'area_code')) %>%
distinct(GSS_CODE, NAME, new_cases, new_cases_per_10k_523_622) #choose which column to be used in the following analysis
#making maps
tmap_mode('view')
## tmap mode set to interactive viewing
breaks_323 = c(10, 15, 20, 25, 30, 35)
tm1 <- tm_shape(Londonborough_merged) +
tm_polygons(
'new_cases_per_10k_population',
breaks = breaks_323,
alpha = 0.5,
palette = 'PuBu',
colorNA = 'white',
title = 'New cases rate'
) + tm_layout(
legend.position = c('left', 'bottom'),
legend.outside = FALSE,
legend.title.size = 1.2,
legend.text.size = 0.75,
frame = FALSE
) + tm_credits('(A) New cases rate in 3.23-4.22', position = c('left', 'top'), size = 1.2) + tm_compass(type = "arrow", position = c("right", "bottom")) +tm_scale_bar(position =c('right', 'bottom'), text.size = 0.6)
#view map 1
tm1
## Credits not supported in view mode.
## Compass not supported in view mode.
## legend.postion is used for plot mode. Use view.legend.position in tm_view to set the legend position in view mode.
#saving map
tmap_save(tm1,'/Users/ziyichen/Desktop/SDSV/GIS/Data_for_CW/Code&data/data/new_cases_rate(323-422).png')
## Map saved to /Users/ziyichen/Desktop/SDSV/GIS/Data_for_CW/Code&data/data/new_cases_rate(323-422).png
## Resolution: 2389.896 by 1845.268 pixels
## Size: 7.966321 by 6.150895 inches (300 dpi)
#making map 2
breaks_523=c(0,0.5,1,1.5,2,2.5,3)
tm2 <- tm_shape(Londonborough_merged_523) +
tm_polygons(
'new_cases_per_10k_523_622',
breaks = breaks_523,
alpha = 0.5,
palette = 'PuBu',
colorNA = 'white',
title = 'New cases rate'
) + tm_layout(
legend.position = c('left', 'bottom'),
legend.outside = FALSE,
legend.title.size = 1.2,
legend.text.size = 0.75,
frame = FALSE
) +
tm_compass(type = "arrow", position = c("right", "bottom")) + tm_credits('(B) New cases rate in 5.23-6.22',
position = c('left', 'top'),
size =
1.2) +
tm_scale_bar(position = c('right', 'bottom'), text.size = 0.6)
tm2
## Credits not supported in view mode.
## Compass not supported in view mode.
## legend.postion is used for plot mode. Use view.legend.position in tm_view to set the legend position in view mode.
tmap_save(tm2,'/Users/ziyichen/Desktop/SDSV/GIS/Data_for_CW/Code&data/data/new_cases_rate(523-622).png')
## Map saved to /Users/ziyichen/Desktop/SDSV/GIS/Data_for_CW/Code&data/data/new_cases_rate(523-622).png
## Resolution: 2389.896 by 1845.268 pixels
## Size: 7.966321 by 6.150895 inches (300 dpi)
#combine two map
t=tmap_arrange(tm1,tm2,nrow=1)
tmap_save(t,'/Users/ziyichen/Desktop/SDSV/GIS/Data_for_CW/Code&data/data/New_cases_rate.png')
## Map saved to /Users/ziyichen/Desktop/SDSV/GIS/Data_for_CW/Code&data/data/New_cases_rate.png
## Resolution: 2100 by 2100 pixels
## Size: 7 by 7 inches (300 dpi)
t
## Credits not supported in view mode.
## Compass not supported in view mode.
## legend.postion is used for plot mode. Use view.legend.position in tm_view to set the legend position in view mode.
## Credits not supported in view mode.
## Compass not supported in view mode.
## legend.postion is used for plot mode. Use view.legend.position in tm_view to set the legend position in view mode.
###Global Moran
#data cleaning for 3.23-4.22 data
london_exclude_city <-
na.omit(Londonborough_merged) #exclude rows contain null data
#data cleaning for 5.23-6.22 data
london_exclude_city_523 <- na.omit(Londonborough_merged_523)
#creating polygon for 323 dataset and 523 dataset
neibour_323 <- poly2nb(london_exclude_city, queen = TRUE)
neibour_323[[1]]
## [1] 19 20 21 22
neibour_523 <- poly2nb(london_exclude_city_523,queen=TRUE)
neibour_523[[1]]
## [1] 19 20 21 22
#assign weight matrix for each neighbouring polygon using row standardisation
lw_323 <- nb2listw(neibour_323, style = "W", zero.policy = TRUE) #each neighbour is assigned a quarter of total weight
lw_523 <- nb2listw(neibour_523,style='W',zero.policy = TRUE)
#computing global Moran's I
moran.test(london_exclude_city$new_cases_per_10k_population, lw_323)
##
## Moran I test under randomisation
##
## data: london_exclude_city$new_cases_per_10k_population
## weights: lw_323
##
## Moran I statistic standard deviate = 3.1736, p-value = 0.0007528
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.35308476 -0.03225806 0.01474307
moran.test(london_exclude_city_523$new_cases_per_10k_523_622,lw_523)
##
## Moran I test under randomisation
##
## data: london_exclude_city_523$new_cases_per_10k_523_622
## weights: lw_523
##
## Moran I statistic standard deviate = 2.5898, p-value = 0.004802
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.28295263 -0.03225806 0.01481388
###Local Moran visualization
#local moran value for data 3.23-4.22
local_moran_323 <- localmoran(london_exclude_city$new_cases_per_10k_population, lw_323)
summary(local_moran_323)
## Ii E.Ii Var.Ii Z.Ii
## Min. :-0.419544 Min. :-0.03226 Min. :0.1102 Min. :-0.71017
## 1st Qu.:-0.002326 1st Qu.:-0.03226 1st Qu.:0.1664 1st Qu.: 0.08039
## Median : 0.149082 Median :-0.03226 Median :0.2155 Median : 0.36073
## Mean : 0.353085 Mean :-0.03226 Mean :0.2460 Mean : 0.82412
## 3rd Qu.: 0.619470 3rd Qu.:-0.03226 3rd Qu.:0.2974 3rd Qu.: 1.59790
## Max. : 2.209853 Max. :-0.03226 Max. :0.4612 Max. : 4.11139
## Pr(z > 0)
## Min. :0.0000197
## 1st Qu.:0.0550665
## Median :0.3595183
## Mean :0.3109317
## 3rd Qu.:0.4679710
## Max. :0.7612010
#local moran value for data 5.23-6.22
local_moran_523 <- localmoran(london_exclude_city_523$new_cases_per_10k_523_622,lw_523)
summary(local_moran_523)
## Ii E.Ii Var.Ii Z.Ii
## Min. :-0.40785 Min. :-0.03226 Min. :0.1106 Min. :-0.80748
## 1st Qu.:-0.05669 1st Qu.:-0.03226 1st Qu.:0.1670 1st Qu.:-0.05456
## Median : 0.16490 Median :-0.03226 Median :0.2164 Median : 0.43173
## Mean : 0.28295 Mean :-0.03226 Mean :0.2470 Mean : 0.68948
## 3rd Qu.: 0.44756 3rd Qu.:-0.03226 3rd Qu.:0.2986 3rd Qu.: 1.17418
## Max. : 1.39046 Max. :-0.03226 Max. :0.4632 Max. : 3.48158
## Pr(z > 0)
## Min. :0.0002492
## 1st Qu.:0.1213157
## Median :0.3343064
## Mean :0.3309740
## 3rd Qu.:0.5217552
## Max. :0.7903052
#plot local moran map
#There are 5 columns of data.
#copy some of the columns (the I score (column 1) and the z-score standard deviation (column 4))
#the z-score (standardised value relating to whether high values or low values are clustering together)
#change local_moran type
local_moran_tibble_323 <- as_tibble(local_moran_323) #change to dataframe
local_moran_tibble_523 <- as_tibble(local_moran_523)
london_exclude_city <-
london_exclude_city %>% mutate(local_moran_I = as.numeric(local_moran_tibble_323$Ii)) %>% mutate(local_moran_z =as.numeric(local_moran_tibble_323$Z.Ii))
london_exclude_city_523 <-
london_exclude_city_523 %>% mutate(local_moran_I = as.numeric(local_moran_tibble_523$Ii)) %>% mutate(local_moran_z =
as.numeric(local_moran_tibble_523$Z.Ii))
#setting color
MoranColours <- rev(brewer.pal(8, "RdBu"))
#plot a map (in Rmd documend, tmap mode changes to viewing, but raw code is plotting mode)
tmap_mode('view')
## tmap mode set to interactive viewing
#plotting local moran map for 3.23-4.22
tm_323 <- tm_shape(london_exclude_city) +
tm_polygons(
'local_moran_I',
alpha = 0.5,
palette = MoranColours,
title = 'Local Moran I',midpoint=NA
) + tm_layout(
legend.position = c('left', 'bottom'),
legend.outside = FALSE,
legend.title.size = 1.2,
legend.text.size = 0.75,
frame = FALSE
) + tm_compass(type = "arrow", position = c("right", "bottom")) +tm_scale_bar(position =
c('right', 'bottom'), text.size = 0.6)+tm_credits('(A) Local Moran in 3.23-4.22',position = c('left','top'),size = 1.2)
tm_323
## Credits not supported in view mode.
## Compass not supported in view mode.
## legend.postion is used for plot mode. Use view.legend.position in tm_view to set the legend position in view mode.
tmap_save(tm_323,'/Users/ziyichen/Desktop/SDSV/GIS/Data_for_CW/Code&data/data/local_moran_323.png')
## Map saved to /Users/ziyichen/Desktop/SDSV/GIS/Data_for_CW/Code&data/data/local_moran_323.png
## Resolution: 2389.896 by 1845.268 pixels
## Size: 7.966321 by 6.150895 inches (300 dpi)
#plotting local moran map for 5.23-6.22 dataset
tm_523 <- tm_shape(london_exclude_city_523) +
tm_polygons(
'local_moran_I',
alpha = 0.5,
palette = MoranColours,
title = 'Local Moran I',midpoint=NA
) + tm_layout(
legend.position = c('left', 'bottom'),
legend.outside = FALSE,
legend.title.size = 1.2,
legend.text.size = 0.75,
frame = FALSE
) +tm_compass(type = "arrow", position = c("right", "bottom")) +tm_scale_bar(position =
c('right', 'bottom'), text.size = 0.6)+tm_credits('(B) Local Moran in 5.23-6.22',position = c('left','top'),size = 1.2)
tm_523
## Credits not supported in view mode.
## Compass not supported in view mode.
## legend.postion is used for plot mode. Use view.legend.position in tm_view to set the legend position in view mode.
tmap_save(tm_523,'/Users/ziyichen/Desktop/SDSV/GIS/Data_for_CW/Code&data/data/local_moran_523.png')
## Map saved to /Users/ziyichen/Desktop/SDSV/GIS/Data_for_CW/Code&data/data/local_moran_523.png
## Resolution: 2389.896 by 1845.268 pixels
## Size: 7.966321 by 6.150895 inches (300 dpi)
#combine two map
t_c=tmap_arrange(tm_323,tm_523,nrow=1)
t_c
## Credits not supported in view mode.
## Compass not supported in view mode.
## legend.postion is used for plot mode. Use view.legend.position in tm_view to set the legend position in view mode.
## Credits not supported in view mode.
## Compass not supported in view mode.
## legend.postion is used for plot mode. Use view.legend.position in tm_view to set the legend position in view mode.
tmap_save(t_c,'/Users/ziyichen/Desktop/SDSV/GIS/Data_for_CW/Code&data/data/local_moran_combined.png')
## Map saved to /Users/ziyichen/Desktop/SDSV/GIS/Data_for_CW/Code&data/data/local_moran_combined.png
## Resolution: 2100 by 2100 pixels
## Size: 7 by 7 inches (300 dpi)
###Geary’s test result
#geary's test
#geary's test for 3.23-4.22 data
geary.test(london_exclude_city$new_cases_per_10k_population, lw_323)
##
## Geary C test under randomisation
##
## data: london_exclude_city$new_cases_per_10k_population
## weights: lw_323
##
## Geary C statistic standard deviate = 2.9614, p-value = 0.001531
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic Expectation Variance
## 0.63209605 1.00000000 0.01543425
#geary's test for 5.23-6.22 data
geary.test(london_exclude_city_523$new_cases_per_10k_523_622,lw_523)
##
## Geary C test under randomisation
##
## data: london_exclude_city_523$new_cases_per_10k_523_622
## weights: lw_523
##
## Geary C statistic standard deviate = 2.1708, p-value = 0.01497
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic Expectation Variance
## 0.73102781 1.00000000 0.01535192
###Getis Ord General G result
#Getis Ord General G for 3.23-4.22
globalG.test(london_exclude_city$new_cases_per_10k_population,lw_323)
## Warning in globalG.test(london_exclude_city$new_cases_per_10k_population, :
## Binary weights recommended (sepecially for distance bands)
##
## Getis-Ord global G statistic
##
## data: london_exclude_city$new_cases_per_10k_population
## weights: lw_323
##
## standard deviate = 2.5827, p-value = 0.004901
## alternative hypothesis: greater
## sample estimates:
## Global G statistic Expectation Variance
## 3.347784e-02 3.225806e-02 2.230541e-07
#Getis Ord General G for 5.23-6.22
globalG.test(london_exclude_city_523$new_cases_per_10k_523_622,lw_523)
## Warning in globalG.test(london_exclude_city_523$new_cases_per_10k_523_622, :
## Binary weights recommended (sepecially for distance bands)
##
## Getis-Ord global G statistic
##
## data: london_exclude_city_523$new_cases_per_10k_523_622
## weights: lw_523
##
## standard deviate = 1.525, p-value = 0.06362
## alternative hypothesis: greater
## sample estimates:
## Global G statistic Expectation Variance
## 3.336104e-02 3.225806e-02 5.230748e-07
###Descriptive statistics summary
library(ggplot2)
#histrogram
#only use rate of new cases to plot histrogram
covid_323_hist <- hist(covid_323_422$new_cases_per_10k_population,col='Dark blue',density=10,xlab='Rate of new cases',main='Rate of new cases during 3.23-4.22',angle=45,ylim=c(0,15))
mean_323 <- mean(covid_323_422$new_cases_per_10k_population)
var_323 <- var(covid_323_422$new_cases_per_10k_population)
std_323 <- sd(covid_323_422$new_cases_per_10k_population)
mean_323
## [1] 23.5875
var_323
## [1] 25.43532
std_323
## [1] 5.043344
covid_523_hist <- hist(covid_523_622$new_cases_per_10k_523_622,col='Dark blue',density=10,xlab='Rate of new cases',main='Rate of new cases during 5.23-6.22',angle=45,ylim=c(0,15))
mean_523 <- mean(covid_523_622$new_cases_per_10k_523_622)
var_523 <- var(covid_523_622$new_cases_per_10k_523_622)
std_523 <- sd(covid_523_622$new_cases_per_10k_523_622)
mean_523
## [1] 1.771875
var_523
## [1] 0.2949899
std_523
## [1] 0.5431297